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Abstract 

We study the equilibrium beiiaviour of a deterministic four-state mutation-selection model 
as a model for the evolution of a population of nucleotide sequences in sequence space. The 
mutation model is the Kimura 3ST mutation scheme, and the selection scheme is assumed to 
be invariant under permutation of sites. Considering the evolution process both forward and 
backward in time, we use the ancestral distribution as the stationary state of the backward 
process to derive an expression for the mutational loss (as the difference between ancestral and 
population mean fitness), and we prove a maximum principle that determines the population 
mean fitness in mutation-selection balance. 

1 Introduction 

The mathematical modelling of populations subjected to the competi ng evolutionary forces of 
mutation and selection has a long and rich history, see, e.g., HienidigTi). The various approaches 
that have been employed to describe DNA evolution at the molecular level can be classified into 
two main categories, comprising stochastic approaches on the one hand and deterministic on the 
other. 

Th e stocha s tic m odels deal with finite populations using Wright-Fisher sampling ijWrightl 
l|l93lj) : llilwenj l)l979j) . Ch. 3). The stochastic formulation lies at the heart of the neutral the- 
ory and had a strong influence o n the methods used today to analyse population sequence data, 
see, e.g., iHartl and ClarkI l)l997(l : iLi and GrauiH |1990) for reviews. Selection can be treated as 
well; however, this is limited to very simple situations with two alleles only. The analysis of more 
complicated settings becomes impractical. 

In the deterministic formulation, more challenging selection schemes can be treated, at the 
cost of neglecting genetic drift. Classical mutation-selection models forrnulate differential equa- 
tions for the evolu t ion of gene frequencies in infinite populations (Fishen l|l922f) : [Haldanc ( 1 928|) : 
ICrow and Kimural l)l970|) . Ch. 6; for an up-to-date review see FBiirge j l)200(]|) ~ These have been 
adapted to sequence space by identifying alleles with s equenc es, and choosing an app ropriate mu- 
tation mechanism — a tradition that was started bv |Eigenj |l971') an d reviewed bylEigen et all 
lll989l). A number of results can be found, e.g., in lO'BrienI l)l985j) : iLeuthausseJ l)l986l Il987|) : 
iRumschitzkl ()l987l) : lTarazonal l)l992(l and, in a time-continuous formulation where mutation and 
selecti on are decoupled from each other, in lBaakd (1995); Baa ke et all l)l997(l : lBaake a nd Wa gneJ 
l)200l)) . But this body of literature almost exclusively deals with two-state models, where each 
site in the sequence can be occupied with one of two states (or alleles), wildtype or mutant, in the 
assumption that this captures the essential behaviour of DNA sequences which are written in a 
four-letter alphabet. 

In a recent article. iHermisson et all l)200 j) refined this approach to a four-state model within a 
physical framework, concentrating on linear and quadratic fitness functions. The results obtained 
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for the four-state model show a much richer behaviour than is observed for the two-state case. 
Therefore, four-state models clearly deserve further investigation. In this article, we consider a 
four-state model for general permutation- invariant fitness functions, formulated entirely within 
the biological framework. Our main result is a maximum principle that allows us to determine 
the population mean fitness in mutation-selection balance by maximising the difference between 
the fitness and a suitably defined function that describes the mutational loss. 

In practice, permutation-invariant fitness means that the fitness only depends on the type and 
the number of mutations relative to a reference sequence, not on their position in the sequence. 
This is, of course, a strong restriction; however, because of the difficult accessibility and complexity 
of re alistic fitness landsca p es, it is a wid ely used assumption in the population genetics literature 
re.g.. lCharlesworthl lll99fl() : IWiehd lll997l) 'l. Furthermore, a permutation-invariant fitness function 
describes the accumulation of many small mutational effects surprisingly well. It is actually a 
good approximation for the fitness functi on in some models for concr e te experimental situation s 
like the DNA-binding models treated by Ivon Hippel and Ber3 l|l98(fl ; 'Gerl and and Hwal l)2002j) , 
which are formulated as two-state models and would certainly profit from a generalisation to a 
four-state description. 

The four-state model of 'Hermiss on et all l|200lt ) is a generalisation of the two-state model or 
biallelic chain ( Baake et ai, 1997). Mutation-selection models of this type are closely connected 
to certain models of statistical physics, so-called quantum spin chains. Whereas the biallelic chain 
is related to the quantum Ising chain, the four-state model considered in this article cor responds 
to the Ashkin- Teller quantum chain, see, e.g., tKoh mo to_et al, (198^ and iBaxterl l)l982|) . Ch. 12. 
However, this correspondence does not mean that results from statistical physics can be transferred 
directly to biology, because some of the quantities considered in t he context of statistical physics 
are not those that are of interest here, compare the discussion bv lBaake and W agner (2001). 

The outline of this paper is as follows. We start with a general introduction to mutation- 
selection models, defining the ancestral distribution and the observable quantities in this class 
of models. This is then specialised to the sequence space of the four-state model which we are 
interested in, choosing appropriate models for mutation and selection. Exploiting the permutation 
invariance of the fitness function and the symmetries of the mutation model, the sequence space 
can be reduced to the permutation-invariant subspace. Subsequently, we define the mutational 
loss function and state the maximum principle. This principle holds exactly in three special cases 
which we discuss in detail. Finally, this is followed by a summary and a brief outlook. The rather 
technical proofs of the maximum principle are given in two appendices. 

2 Mutation-selection models 

We consider a population of haploid individuals^ whose genotypes are chosen from a sequence 
space, a set of a finite number i' of possible genotypes i. The population is described by the 
population distribution p, a i/-dimensional vector with entries Pi > 0, indicating the relative 
frequency of type i in the population. Hence, p has to be normalised such that X^iLiPi ~ ^■ 
For finite population size, the pi are rational numbers. In this article, we concentrate on the 
deterministic limit of an infinite population size, where the relative frequencies can take real 
values. 

Ignoring environmental effects, mutation and selection are assumed to depend only on the 
genotypes of individuals. In this framework, the evolutionary processes to be considered are birth 
and death of individuals, and mutation from one type to another. 

In the time-continuous model, each individual of type i gives birth to an identical copy with 
a rate bi and dies with a rate d,;, hence we have an effective reproduction rate ~ bi — di, also 
called the Malthusian fitness (Biirger (200(1), Ch. 1). These values are collected in a diagonal 
reproduction matrix R — diag(ri, . . . , r^). 

^The theory applies as well to populations of diploids without dominance, where the evolution equations reduce 
to those of a haploid population (cf. iBiirger. 1,2000.) . Ch. 2.2) 
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Mutation from type j to i occurs with a rate My . To preserve the normahsation of the 
population distribution p, the diagonal entries Ma of the mutation matrix M = (Mij) are chosen 
such that M has a vanishing sum over the columns, J^i ^ij — 0, which makes M a Markov 
generator. Unless we talk about unidirectional mutation (cf. Sec. 15. we will assume that the 
mutation matrix is irreducible, i.e., each genotype can be reached from any other by mutation, 
possibly in several steps. With the definition of the time-evolution operator H = R + M, this 
leads to the evolution equation 

p{t) = {H-r{t)\)p{t), (1) 

w here r is the population m ean fitness r (t) — Pi(t) and 1 denotes the v y. v identity matrix, 
cf. lCrow and Kimural lll97f]l) . Ch. 6, and lBikg er (2000), Ch. 3. 

Irreducibility of M implies that of H , and the Perron- Frobenius (PF) theorem guarantees that 
there exists a unique stable equilibrium solution, which is given by the strictly positive eigenvector 
p corresponding to the largest eigenvalue Amax of H , 

Hp = AniaxP- (2) 

In the limit as t — > oo, the population distribution converges towards this equilibrium solution, 
linit^ooP(0 -■ P- 

2.1 The ancestral distribution 

We are particularly interested in the equilibrium solutions p — Q. In this case, Eq. becomes an 
eigenvalue equation for H with eigenvalue Amax = r. The right PF eigenvector p is the population 
distribution in equilibrium, whereas the entries Zi of the le ft PF eigenvector 2 : deter mine the relative 
reproductive success of type-* i ndividuals, as sh o wn by iHermisson et al\ l)2002() . The ancestral 
distribution, also introduced bv IHermisson et "all l)2002j) . is a probability distribution defined as 
ai = ZiPi, with the normalisation of z chosen such that = 1. Here, specifies the fraction 

of the equilibrium population whose ancestors, an infinitely long time ago, were of type i. 

In analogy to the way that the population distribution is defined as a time-dependent quantity, 
this can also be done for the relative reproductive success z and the ancestral distribution a. 
However, this demands some notational efforts, and as we do not need this property later on, we 
limit ourselves to the definition of the ancestral distribution as an equilibrium quantity. 

2.2 Means 

A population is macroscopically described by mean quantities. We introduced two probability 
distributions, hence there are two types of averages that are relevant in our model. Every mapping 
that assigns a value Oi to each possible genotype i can be averaged with respect to the population 
distribution or the ancestral distribution. 

The population mean of o, denoted by o(t) , is given by 

o(t) := J2o,p,{t). (3) 

i 

The population mean in equilibrium, i.e., in the limit as t — > 00, is denoted by o. 
The ancestral mean of an operator o, denoted by o, is defined as 

o ^OjCi . (4) 

i 

Note that the ancestral mean does not depend on time, as we defined the ancestral distribution 
as an equilibrium quantity only. 
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3 The four-state model 



The genetic information is coded in the DNA as a string composed of the purines adenine and 
guanine {A,G) and the pyrimidines cytosine and thymine {C,T). 

We consider DNA strands of f ixed l ength TV, which may, for instance, code for an enzyme, 
as modelled by iHermisson et "all (l200lh. The f our b asic states {A, G, C, T} are mapped onto 

{0,1,2,3}, or, as it is done bv IHermisson et all ()200l|) . onto {++, H — +, }■ Conveniently, 

one can exploit the freedom in the choice of this mapping by introducing a relative rather than 
an absolute mapping between the bases {A, G, C, T} and the symbols {0, 1, 2, 3}. This essentially 
means that one can choose independent mappings at each position along the DNA strand. At 
any position, the mapping can be defined such that the symbol 0, or ++, corresponds to the 
corresponding base in a given preferred sequence, which is usually chosen to be the wildtype or 
master sequence of ma ximal fitness rma v - Thu s the wildtype sequence is mapped onto the sequence 
{0}^ = 000... 0, see IHermisson" et cfif 1)200 ill for details. The mapping between the remaining 



nucleotides and the symbols 1,2,3 will be discussed below. The sequence space consists of all 
possible A^-letter sequences in these four symbols, so it is given by {0, 1, 2, 3}^ and has dimension 

In what follows, we shall not really need the complete information about the sequences. It will 
be sufficient to characterise a sequence by its mutational distance with respect to the wildtype 
sequence {0}^, which just counts the deviations from the wildtype sequence. Whereas in the 
two-state model the mutational distance is given by a single integer, which counts the number of 
bases along the DNA strand that differ from those in the wildtype, we now need three non-negative 
integers di, d2 and da, according to the different types of mutations that can occur. We define the 
mutational distance d of a sequence as 

d = id,] := [ #(2) , (5) 
\d3j V#(3)/ 

where #(1), #(2) and #(3) denote the number of entries 1, 2 and 3 in the sequence, respectively. 
The total mutational distance is defined as the sum d := di+d2+d3, which takes values < d < iV. 



3.1 Mutation 

Mutation is taken to be a point process that acts at each site independently. Disregarding more 
complicated mechanisms such as deletions and insertions, we only take into account the replace- 
ment of one base by another. Taken over the whole sequence, this happens with certain rates fj,k, 
where k indicates the type of replacement. We allow only one mutation at a time, as modelled 
by a P oisson process. This leads to a single step mutation model, which was first introduced 
bv lOht a and Kimura ( 197 3J). We wor k with the Kimura 3ST m utation scheme shown in Fig. ^ 
l)Kimural l|l984 : ISwofford e.t al\ jl99fil) : lEwens and Cxrantl l|200lh . Ch. 13) which assumes that, of 
a possible 12 mutation rates that can be chosen in our setting, only three different mutation rates 
Hi, 1^2 and /X3 occur. In particular, forward and backward mutation rates are the same, and the 
mutation process respects a symmetry between exchanges of purines and pyrimidines. 

This mutation scheme can be treated to various degrees of sophistication. Apart from the 
full Kimura 3ST scheme, where all three mutation rates are different, there are also two simpler 
models that are worth mentioning. 

The simplest approach is to take all mutatio n rates to be equal, /x i ^ fJ.2 ~ fJ-a- This case is 
known as the Jukes-Cantor mutation scheme (Jukes and Canton Il969|) . 

Due to the similar shapes of the nucleotides, transitions, i.e., the replacement of one purin/ 
pyrimidine by the other, are more frequent than transversions, i.e., the replacement of a purin/ 
pyrimidine by a pyrimidine/purin. The mutation rates describing the transversions are fairly 
similar, so /ii « /X3, whereas the mutation rate for the transitions /i2 is typicall y larger by a fac tor 
of about 2 to 40. This is taken into account in the Kimura 2 parameter model l|KimuraL ll98Cl|) by 
assuming that jJ-i — ^3 < ^i2- 
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Figure 1: The Kimura 3 ST mutation scheme. 



Assume there is one particular sequence sq with maximal fitness rmax, the wildtype or master 
sequence, which is mapped onto {0}^. For any other sequence, the corresponding representation 
in terms of the symbols {0, 1, 2, 3} is then obtained by comparing it to the wildtype, and assigning 
one of the labels 1, 2, 3 at each position where it differs from the wildtype sequence, according to 
the type of mutation as given in Fig. ^ Analogously to the mutational distance, we can define 
the Hamming distance |Hammiiig (1950); van Lint (1 982) , Ch. 3) between two sequences Si and 
Sj, by comparing the sequences with each other. The restricted Hamming distances dk{si, Sj) are 
the numbers of type-fc mutations between the sequences Si and Sj, i.e., mutations that occur with 
rate ^k] the total Hamming distance is d{si, Sj) = (ii(si, Sj) + d2{si, Sj) + d3{si, Sj). 

In the Kimura 3ST setting, the entries Mij of the mutation matrix are given by 



Mfc 
N 




for d{si, Sj) = dk{si,Sj) = 1 , 
for d{si, Sj) > 1 , 

3 

^ fik for i = j , 



(6) 



fe=i 



where, as mentioned above, the diagonal entries are chosen such that M is a Markov generator 

Here, the mutation rates are scaled as mutation rates p 

whole DNA string being constant, see the discussion in lBaake and Wagned l)200ll 



er sitej_with the mutation rate over the 



3.2 Selection 

Whereas the process of mutation is well understood and straightforward to model, the choice of 
the fitness landscape on the molecular level is far from being clear. Realistic fitness landscapes 
would be rather rugged and strongly dependent on the function of the DNA sequence, but they 
are hard to access experimentally. 

We shall use the severe simplification of a permutation-invariant fitness function, which is nev- 
ertheless a rather common (and usually implicitly made) assumption in theoretical investigations 
of mutation-selection models such as Charlesworth (1990); Wiehc ( 19 9^ as well as in the mod- 
elling of concrete e xperim ental settings like in DNA-binding models (jvon Hippel and Berel Il986t 
iGerland and Hwal . l2002() . Using permutation- invariant fitness, one assumes that the fitness of a 
sequence depends only on the number of mutations of the various kinds, not on their location 
within the sequence. Hence, we can describe a sequence completely by its mutational distance 
d = (^1,^2,^3) with respect to the wildtype sequence. As there are 4^ different sequences, but 
only {N -f 1)(7V -I- 2){N + 3)/6 different distances d with < d < N, this reduces the effec- 
tive type space enormously. In the permutation-invariant fitness model, the number of possible 
different genotypes, and thus the dimension of the permutation-invariant subspace, is given by 
1/= {N +1){N + 2){N + 3)/6. 
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4 Reduction to the permutation-invariant subspace 



In our model, three different spaces are relevant: (i) the 4^-dimensional full sequence space, (ii) 
the reduced sequence space of dimension v, and (iii) the three-dimensional space of the mutational 
distances. 

In the full sequence space of dimension 4^, each sequence corresponds to a different basis 
vector, and the population p is then completely determined as a point on the (4^ — l)-dimensional 

hyperplane defined by X]i=i?'i ~ where the projection on each axis gives the frequency pi of 
the corresponding sequence. 

Analogously, the j/-dimensional reduced sequence space, which is the permutation-invariant 
subspace of the full sequence space, is spanned by unit vectors, each of which corresponds to 
the set of all sequences which have the same number of mutations of each type, i.e., to one of 
the u different mutational distances d. Here, the population p is also given as a point on a 
hyperplane (of dimension v — 1). In general, the transition from the full to the reduced sequence 
space is accompanied by a loss of information. As long as we consider systems with a unique 
equilibrium population, we know that this equilibrium will be permutation-invariant, because 
starting from a permutation-invariant initial population, we will reach a permutation-invariant 
equilibrium because the fitness function, as well as the mutation scheme, disregard the order in the 
sequences. As the equilibrium is unique, it will be reached from any initial condition. Therefore, 
sequences with the same numbers of each type of mutation, i.e., the same d, must occur with 
the same frequency in the equilibrium population. Thus, as long as we are only interested in 
equilibrium properties of systems with a unique equilibrium, it suffices to restrict ourselves to the 
reduced sequence space. 

Finally, we have the three-dimension mutational distance space of the mutational distance 
vectors d with Cartesian coordinates di, d2 and d^. The basis of this space is formed by the 
Cartesian unit vectors ei = (1,0,0)*, 62 = (0,1,0)* and 63 = (0,0,1)*, which are the basic 
directions of mutation. The condition Q < d — di -\- d2 -\- d^ < N restricts the possible mutational 
distance vectors d to a simplex in the positive quadrant, as shown in Fig. |21 There is a one- 
to-one correspondence between the sequences in the reduced sequence space and the mutational 
distance vectors d in this simplex, and we label the elements of the reduced sequence space by the 
corresponding mutational distance vectors d. Whenever we speak of a sequence d, we refer to the 
corresponding sequences in the reduced sequence space. 




Figure 2: Mutational distance space in the case of permutation-invariant fitness. 

The mutational distance space is required to define the neighbourhood of sequences. As we use 
a single step mutation model, each sequence, now labelled by d, has at most 12 neighbours, i.e., 
sequences to which they can mutate within a single mutation step. The neighbours of a sequence 
d have mutational distances d ± e^, where ^ determines the direction of mutation and the are 
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combinations of the basis vectors e/j. More precisely, there are two types of mutational directions 
^; firstly, the mutations from wildtype to mutant, ^ = k with k e {1,2,3}, which correspond to 
unit vectors = e^, and vice versa, which correspond to = — efe. For the second type, where 
one type of mutation is replaced by another mutation, one mutation step corresponds to a vector 
= efe — ee. These directions are labelled by pairs ^ = {+k, —i) with k,£ £ {1, 2, 3} and k > i, 
and the corresponding mutations in the inverse directions are labelled accordingly with k < i. 
Finally, we note that points on the surface of the simplex in mutational distance space have fewer 
neighbours. Clearly, only those mutations that do not leave the simplex are permitted. 



4.1 Similarity transformation 



We exploit the permutation invariance of the fitness function to reduce the sequence space to its 
relevant permutation-invariant subspacc. Therefore, wc transform our time-evolution operator H, 
which is a 4^ X 4^-matrix, to a matrix Hpi^ of dimension = {N + 1){N + 2){N + 3)/6 that 
describes the reduced sequence space only, the subscript "piv" refers to the permutation-invariant 
subspace. This is done by the means of a similarity transformation T. We have 



T-^HT 



H 



piv 





(*) 
H' 



(7) 



where ffpiv is the i^-dimcnsional time-evolution operator describing the reduced sequence space. 

A condition on the transformation T is that it preserves the Markov property of the mutation 
matrix M. Hence, it must be an L^-transformation, which is guaranteed by the property = 
1. For the relevant permutation-invariant subspace, wc need to combine all sequences that belong 
to the same mutational distance vector d. This, together with the Markov condition, determines 
the entries in the first u columns of T. In the column assigned to d, each entry is either 0, or, if 
the sequence corresponding to the row in question has the mutational distance d, it is l/n^, where 
na is the number of sequences that are mapped onto d. This number is given by the multinomial 
coefficients 

nd = \ , , , J I = ,,,,,,,, , (8) 



do,rfi,rf2,d3 



do\di\d2\d3\ 



with do := X)fe=i dk denoting the number of wildtype sites. The actual choice of other columns 
of T does not influence the submatrix iifpiv, here we only need that T is invertible. 

For a sequence length of N = 2, for example, the relevant part of the transformation T has 
the form 



000 100 010 001 200 110 101 020 Oil 002 
/ 1 
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1/2 
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1/2 
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1/2 
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1/2 



1/2 
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33 



(9) 
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where the triples at the top give the mutational distances d to which each column corresponds, 
whereas on the right the actual sequences s corresponding to each line are displayed. Only non- 
zero entries are shown. In this case, the remaining six columns of T, shown symbolically as (*), 
correspond to the antisymmetric subspace, which contains sequences with mutational distances 
d = (100)*, (010)*, (001)*, (110)*, (101)*, (Oil)*. 

The diagonal entries of .ffpiv remain unchanged compared to the original H; they are Hpi^M = 
ra — X^feMfc- The off-diagonal entries, i.e., the mutation rates u in the permutation-invariant 
subspace, depend on the direction of mutation. Using the normalised versions of the mutational 
distances Xk ■= dk/N, they are given by 

d^ d + Bk ■■ u^^ = jjLk xo (3 eqns.) 

d^d-ek '. Wrf = fikXk (3 eqns.) (10) 

d^ d + ek- ee : u^'''~^ = l^m xe (6 eqns.) 

where k,l,m G {1,2,3} are pairwise different, so {k,l,m} = {1,2,3}. Our notation is such 
that u'^'' and u^*^ denote the rates for mutations from distance d in the positive and negative 
k direction, respectively, and u^''~^ the corresponding rate in direction Cfe — e^. The mutation 
rates now depend on d, reflecting the fraction of sites that can mutate with the specifled effect. 
In particular, note that this implies that iJpiv is not a symmetric matrix. The above definition of 
mutation rates takes care of the boundary condition u^^ = for d on the boundary of the relevant 
simplex in the mutational distance space with ±ej pointing outwards. 



4.2 Symmetrisation of the time-evolution operator iifpiv 

There is an alternative way to arrive at a matrix that describes the permutation-invariant subspace, 
which, however, does not preserve the Markov property for M . As the original 4^ x 4^ matrix 
H is real symmetric, we can block-diagonalise it by an orthogonal transformation O. This is an 
L^-transformation, and it preserves the symmetry of the matrix, so the corresponding vxv matrix 
i?piv given by ^ 

O'HO = (^^^^ 1^ (11) 

is symmetric, where O* = denotes the transpose of the orthogonal matrix O. In this case, 
the states in the permutation-invariant subspace are again superpositions of all sequences of equal 
distance d, but now with coefficients 1/^/nd, as opposed to l/rid for the -transformation T. 
Knowing this connection, we can symmetrise the permutation-invariant part iipiv of the time- 
evolution operator. By slight abuse of notation, we get for the permutation-invariant part 

Hpi. = (O'HO)^.^ = (O'TT-'HTT-'O)^.^ = D^Hpi^D, (12) 

where D := (T^^O)piv, and the subscripts mean that we restrict to the permutation-invariant 
part. As the relevant columns of the matrices T and O differ only by factors ^/nd, the trans- 
formation D that symmetrises the time-evolution operator iJpiv is in fact diagonal, with entries 

Ddd = \frid- 

The corresponding mutation rates u of the symmetrised system are given by 





1 nd 






= 


/ nd 


"d 


V ^cl-e,. 


■+k-t _ 




U 



lik\jxo{xk + ^) = \fu^^u^l^^, 



rid+e 



= Hm\lxe{Xk + —) = \Ju^ Md+e,-e, • K^^) 



8 



The property that Hpi^ is symmetrisable by means of a diagonal transformation allows us 
to write the eigenvalue equation ^ in ancestral formulation, which is the starting point for the 
proofs of the maximum principle. In fact, the proofs presented below can be formulated for any 
mutation matrix that is symmetrisable by a diagonal transformation, which is equivalent to the 
property that it describes a reversible process, 

M.jQj = Mj,q,, (14) 

where q is the equilibrium distribution of the mutation process without selection. In this case, 
we can use the diagonal transformation Qij = Sij^^, and from H14|l . we obtain the symmetrised 

mutation rates My = My \/qjJqi = MijMji. Therefore, the mutational distances are not bound 
to be the numbers of mutations of a DNA sequence, but the model can be reinterpreted in the 
multilocus model context, where the genetic distances are three arbitrary traits that determine 
the fitness. In fact, the maximum principle discussed below can also be derived for a model with n 
states at each site of the sequence, which then has an interpretation as n traits contributing to the 
fitness, as long as we talk about a n-dimensional single step model with a mutation matrix that 
describes a reversible process. A more general and more systematic approach thus seems feasible, 
it will be described bv lBaake et al. (2003). 

In what follows, we drop the subscript, and use H and H to denote the submatrices corre- 
sponding to the permutation-invariant subspace. Both time-evolution operators H and H have 
the same eigenvalues, but the eigenvectors differ. The left and right eigenvectors 5 and p oi H, 
for the largest eigenvalue, are related to the corresponding eigenvectors z and p of iJ by 

z = zD and p = D^p . (15) 

The relation between the ancestral distribution a and the symmetrised population p is given by 

ai = ZiPi = (zD~^)^{Dp)^ = %pi r-^ pf , (16) 

as 2 ~ p due to the symmetry of H . With the relation p ~ -^a, the eigenvalue equation of H in 
ancestral formulation becomes H^/a — f^/a^ which explicitly reads 

for some distance d. Here, ^ determines the six possible directions of mutation, and the sign 
indicates the forward and backward direction. 



5 Maximum principle 

For the permutation-invariant system, we can derive a maximum principle for the population mean 
fitness that involves maximisation only over the three components oi x = djN ai the mutational 
distance space, as opposed to the maximisation over the i^-dimensional reduced sequence space 
according to Ray l eigh's principle. It finds its analogue in the scalar maximum principle given in 
l|Hermisson et ad . l2002l Eqs. (30) and (33)). In the four-state model, we have 

r = sup {r(x) — g{x)) (18) 

X 

with the mutational loss function g{x) defined as 

gix) := J2 ^w+«(a;) + ^-^(a;) - 2^Ju+i{x)u-^x) ^ (19) 
with summation over all six directions of mutation ^. 
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This is exact for three special cases, naraely (i) for unidirectional mutation, (ii) for linear fitness 
and mutation functions, and, most importantly, (iii) in the limit of infinite sequence length. These 
cases will be explained in what follows, and the derivation of the maximum principle for each case 
is presented as well. For other systems, the maximum principle gives an approximation, which, 
under reasonable assumptions, one might expect to differ from the true result by correction terms 
of order 

If the supremum in Eq. H18|l is assumed at a unique value x, which is the generic case, this 
value is the ancestral mean mutational distance ir, and we have, in addition to Eq. H18|l . 

r = r{x) — g{x) = r — g{x) , (20) 

which again is exact for the three special cases mentioned before. This maximum principle, in the 
terminology of physics, is akin to the principle of minimal free energy. 

5.1 Unidirectional mutation 

We now consider the first of the three situations where the maximum principle is exact, the case 
of unidirectional mutation. 

Unidirectional mutation means that only such mutations that increase d happen, the mutation 
rates towards types with smaller or equal d are zero. In the most important case of a monotonically 
decreasing fitness function, this means that all mutations are deleterious. This mutation scheme 
does not go in line with the mutation rates for the DNA system as given in Eq. H10|) . but it is 
equivalent with a collapse of the mutation scheme as shown in Fig. |21 




Figure 3: Simplified mutation scheme for unidirectional mutation. 

For unidirectional mutation, the mutation matrix M is no longer irreducible. Hence, in this 
case, the Perron-Frobenius theorem does not apply. The equilibrium is not unique, but depends 
on initial conditions. Once the wildtype is lost in the population, it can never occur again, because 
the mutation rates back to types with smaller d are zero. The structure of the mutational distance 
space is shown in Fig. 01 where the wildtype on the top corner of the mutational distance space 
"feeds" all mutants underneath. 

Although unidirectional mutation rates do not represent the mutation model for DNA se- 
quences as set up in this article, they still are a reasonable approximation. If in the DNA model 
the selection is sufficiently strong, most individuals present in the population will have a geno- 
type with a small mutational distance from the wildtype. The mutations that leave d constant 
(e^ — ei — Ek), and those that decrease d (e^ = ~ek), which in the case of a decreasing fitness 
function are neutral and advantageous, respectively, occur with rates proportional to Xk, and 
therefore are small for individuals with small mutational distances x = X]fc=i = d/N, which 
form the main part of the population, whereas the mutations that increase d (e^ = e^) happen 
with rates proportional to 1 — x, which are of order 1 for small x. Therefore, it is reasonable to 
approximate the mutation rates of the DNA mod el b y unidi r ection al mutation rates, which is the 
well known infinite sites limit fsee iKimural l)l969^ or lEwenj l)l979j) . Ch. 8). 

In the case of unidirectional mutation, the mutational distance space can be divided into three 
domains with respect to each sequence d, namely the ancestral cone, the offspring cone and the 
sibling domain. Here, all sequences that can mutate to d lie in the ancestral cone AC{d), all 
sequences that d can mutate to lie in the offspring cone OC{d), whereas the sequences that are 
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Figure 4: Structure of the mutational distance space for unidirectional mutation. 



not connected to d via a mutational path form the sibling domain SD{d), which is the remainder 
of the mutational distance space. 

The eigenvalue equation of H for unidirectional mutation is given by 

(f-Xd)Pd = Xl^^-e,Pd-e, with Ad = ra-gixd) = - (^1) 

k k 

where the diagonal entries are the eigenvalues Ad, because H is a lower triangular matrix. In 
this case, the situation is particularly simple, and we can directly infer some properties of the 
population distribution p, which allow us to show that the maximum principle holds in this case. 
Suppose there is one d such that p-^ > 0, but Pd-e^ = for /c e {1, 2, 3}. This can only happen 

if f = A^. An evaluation of Eq. (|21(l for d — Sk and for other sequences in the ancestral cone 
yields that no sequence in the ancestral cone can contribute to the population, i.e., Pd = for all 
d e AC(d). 

If there was a sequence d^ in the offspring cone of d with Ad+ > A^, we would get Pd+ < Oj 
which is a contradiction to the condition pi > 0. Thus, Ad < A^ for all d E OC{d), which 
corresponds to the maximum principle l|18() as r = max^^^p^^^ [r{xd) — gixd)]. If now A^+ = 

A^ — "^"3 get Pd — for all d e AC{d'^), including d, so that in this case the offspring cone 
OC(d^) spans the population rather than OC{d). Evaluating the eigenvalue equation for the 
sequences in the offspring cone, we get Pd > for all d e OC{d) or d S OCld^), respectively. 

All sequences in the sibling domain SD{d) descend originally from the ancestors of d which 
are not present in the population. Thus, their frequencies must vanish, unless there is a sequence 
d^ e SD{d) with Ad+ ~ Aj. In this case, if Pd+ > 0, the sequences in both offspring cones have 
non- vanishing frequency pd > for all d G OC{d)UOC{d'^), the frequencies of all other sequences 
vanish. 

As the mutation matrix is not irreducible for unidirectional mutation, the population distribu- 
tion in equilibrium is not unique, but depends on the initial conditions. The equilibrium with the 
highest mean fitness is always assumed if the wildtype initially occurs with non-zero frequency, 
or if we consider the limit of small, but non-vanishing back mutations. In this case, we have 
r = Amax = supd [rixd) - g{xd)]- 

If, however, the wildtype is not present in the initial population p{t = 0) and the mutation 
rates for constant or decreasing d are exactly zero, we have to consider only the part of the 
mutational distance space that is spanned by the offspring cones IJ^. OC(din) of all initially 
present sequences din. Now, the mean fitness assumes the highest possible value in this subspace 
f = supy Q(-;(^. J [i~{xd) ~ 9ixd)\, which is assumed for at least one sequence . 
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If this maximum is unique, the equihbrium population is given by the right eigenvector corre- 
sponding to this eigenvalue , which has non-zero entries only for the sequences in the offspring 
cone OC(<i^). The left eigenvector corresponding to has non-zero entries only for the ances- 
tors of , so that the only d with Od 7^ is d^, and hence = d is the only ancestor. This 
yields Eq. (j23). 

An interesting case arises, however, if the maximum is not unique, but attained at two sequences 
d and d^ in the subspace under consideration. In this degenerate case, the ancestral distribution 
cannot be obtained as easily as shown above, as the left and right eigenvectors have no non-zero 
overlap and thus it is not possible to normalise z such that Ui = 1. 

We now have to distinguish between two cases. If the sequences lie in parent-offspring relation, 
i.e., d^ e OC{d), d is the single ancestor, whereas the population is formed by the offspring cone 
of d^, i.e., > if and only if d e OC(d"*"). Note that Eq. if^ still holds, although in this 
special case the only ancestor d has zero frequency in the population. 

If, however, d and d^ lie in sibling relation to each other, the population is formed by the 
unification of their offspring cones OC(d) U OC(d+); and d and d+ both may have non- vanishing 
ancestral frequency, which are then determined by the initial conditions. 



5.2 Linear model 

We now consider the second case where the maximum principle applies exactly. Here, both the 
fitness function and the mutation rates depend linearly on functions of the genotype components 
X]~, and thus can be written as 

r{x) = ro-^akyk(.Xk) and w^^{x) = + ^(3^^ yk{xk) , (22) 

k k 

with p arar neters ak and ■ This type of model has been used, e.g., in Ivon Hippel and Berd 
l|l98fil) and lCerland and Hwal 1 20021) in a two-state version. 

In this case, the mean fitness can be obtained by maximisation over the three components of 
a;, and the maximum principle in the form of Eqs. (|18|l and H20I) holds true for yk{xk) = Xk- This 
can be shown by a direct calculation starting from Eq. (|17|l . which is carried out in Appendix 1X1 



5.3 Infinite sequence length 

Finally, we consider the limit as the sequence length becomes infinite. In the limit as TV — > cxd, we 
use Xd — d/N to describe the mutational distance to the wildtype. The quantities Xd^ — dk/N 
and Xd = X]fc=i fulfill the inequalities < Xd,, < 1 and < a;^ < 1- As long as we operate 
with finite N, the Xd^. take discrete rational values only; in the limit N — > 00, however, they 
become dense, and it is thus reasonable to pass to a continuum formulation. Consider the fitness 
function r{x) and the mutation rates u^^{x) as functions defined on the mutational distance 
space. The mutation rates are positive, continuous functions, obeying the boundary condition 
that they vanish for all x at the boundary of the simplex in the mutational distance space, where 
they correspond to mutations out of the set of possible mutational distances. The fitness function 
has to be piecewise continuous, i.e., it can have discontinuities only along (finitely many) surfaces 
in the mutational distance space; at the discontinuities, it must be either left or right continuous. 
This allows for a number of biologically meaningful fitness functions, like, for example, truncation 
selection. For any finite N, the fitness and mutation functions are sampled at all possible Xd- The 
limit iV ^ CX3 is carried out such that the functions are kept constant, but the sampling gets finer 
with increasing N. 

In the limit as ^ 00, Eqs. (|18(l and H20|l hold true. The proof is presented in Appendix IbI 
Even for rather short sequence lengths, the maximum principle yields a reasonable approximation. 

Figures [S] and show numerical results obtained using the maximum principle in comparison 
with those for finite systems obtained by a direct diagonalisation of the time-evolution operator. 
For both figures, a quadratic symmetric fitness function r{x) = (1 — X] ^fe)^ been used, but they 
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Figure 5: Population mean fitness r and ancestral mean mutational distance Xk for vary- 
ing mutation rate fi for a model with Jukes-Cantor mutation scheme and quadratic fitness 
r = (1 — The curves for finite sequence lengths N = 5,10,20 are obtained by di- 

rect diagonalisation of the time-evolution operator, the curve for A'^ = oo is the result of the 
maximum principle. 



differ in the mutation scheme. The mutation scheme used in Fig. [Slis the Jukes-Cantor mutation 
scheme, whereas in Fig. El we used a Kimura 2 parameter (K2P) model with /ii = /ia := /i = jqI^2- 
For the totally symmetric Jukes-Cantor model, the ancestral means of the mutational distances 
all coincide, whereas they differ in the K2P model, according to xi — <X2- 

6 Summary 

In this article, we investigated the mutation-selection balance in the mutation-selection model 
introduced by Hcr misson et ai (2001). There, a deterministic approach to model the DNA evo- 
lution of asexual populations was taken. We consider four-state sequences subject to the forces of 
mutation and selection. For simplicity, selection is taken to be permutation invariant, which leads 
to a three-dimensionally structured mutational distance space, and the mutation model is a single 
step model on this structure. 

Using the concept of the ancestral distribution, as introduced bv^Herm isson et~al\ 1)20021) . we 
derived a maximum principle for the population mean fitness r and the ancestral mean genotype x 
in equilibrium, which involves a maximisation over the three dimensions of the mutational distance 
space. 

This maximum principle gives the exact mean fitness in the three limiting cases of (i) unidirec- 
tional mutation, (ii) linear fitness and mutation functions, and (iii) the limit of infinite sequence 
length. For finite sequence lengths, it is an approximation which we expect to be correct up to 
terms of the order 1/A^. Numerically, we found that already rather small sequence lengths were 
well reproduced. 

The maximum principle generalises the results of lHermisson et all l)2002|) . where the case of a 
one-dimensional mutational distance space was treated, which can be interpreted on the level of 
DNA sequences as a two-state model, with states representing wildtype and mutant. In that case, 
r, X and x = r~^(r) could be obtained by a maximisation over one dimension. In our model, we 
have to maximise over the three dimensions of the mutational distance space to obtain r and x, 
whereas the population mean genotype x cannot be derived as easily, because the fitness function 
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Figure 6: The quantities r and Xk for varying mutation rate jjL := fii = /X3 for a model with K2P 
mutation scheme, with ^2 — 10^, and quadratic fitness r = (1 — X]fc=i ^^fc)^- The curves for finite 
sequence lengths = 2, 4, 8 are obtained by direct diagonalisation of the time-evolution operator, 
the curve for N — 00 i& the result of the maximum principle. 



is not uniquely invertible in three dimensions. 

Other quantities o f interest are the c orresp onding variances. The expressions for the variance 
of the fitness given by iHermisson et al\ l|2002fl can only be generalised to our model in the linear 
case, not for the case of infinite sequence length. Neither can the variance of mutational distance 
be obtained in a simply way, because this involves inversion of the fitness function, which does not 
have a unique solution in more than one dimension. 

Although our setup is motivated by a model for DNA evolution, it is valid for a system where 
the fitness depends on three arbitrary traits di, d2 and da with a single-step mutation model, as 
long as the mutation matrix describes a reversible process, i.e., MijQj = MjiQi with q being the 
equilibrium distribution of the mutation process without selection. In fact, the maximum principle 
can be generalised to an n-state model, where we have n—1 different traits determining the fitness. 
In this case, we have to maximise over these n—1 quantities. 

Similarly, the restriction to permutation-invariant fitness functions could be dropped. In this 
case, however, we would have to maximise over the N sites of the sequence, so that it loses its 
use, which lies primarily in the simplicity. 



A Proof of the maximum principle for the Unear model 

Starting with the eigenvalue equation of the symmetrised time-evolution operator H (|17|l 



k,l 
k^l 

we make an ansatz for the Ud such that 

= 

ad 



-k 



(23) 



(24) 
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This is equivalent to 



ad 



1 u 



+k 



and 



dd 



Ck u 

Ci u 



-k+l 
d 
+ k-l 
d—Bk+ei 



(25) 



The latter can be seen using the condition for reversibility (|14|) . To determine the constants Ck, 
we multiply Eq. H24|l by its denominators and sum over all d, which yields the ancestral means of 
the mutation rates 



(26) 



Now, we divide Eq. (^5)1 by y/od: and insert the ansatz and (|^ . 



= rd - ("d 



k+l 



(27) 



Multiplication by ad and summation over all d yields, using the explicit form H2()|) of the Ck, 



r = r 



- ^(m+^ + u^^) + 2 ^ ^u+^u-k + E^ 



k,l 



(28) 



So far, we did not use linearity. If r and u depend linearly on some functions yk{xk), we have 
r — r{y) and = u^^{y). With the definition of the mutational loss function 119f) and in the 
case of yk{xk) — Xk, this is Eq. 

In order to obtain the supremum condition H18() , we consider Eq. (|27() for two different sequences 
d and d' and take the difference, using the explicit representation of fitness and mutation functions 
given in Eq. 



« ^ S - 



E 



-k+l 



k 



■13+^ 

■j+k ' rn 



(29) 



This is just the condition 



= Efl — ['"(y) "5(y)]^^s (2/m(a:m) -2/™(a;;„)) 

which has to be fulfilled for arbitrary x and x' . Hence, we have 

d 



= 



[r{y)-9{y)\y^y for me {1,2, 3}. 



(30) 



(31) 



This is a necessary condition for the existence of an extrcmum at y. A sufficient condition for the 
existence of a maximum of the function r — g in y is that the Hessian 



d^jrjy) ~ gjy)) 



(32) 



-I v=v 
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of the second derivatives in the point y is a negative definite matrix. We have 

n[x) = - Y,cAx)U,{x)Ul{x), (33) 

z 

with c,{x) = \{u+^{x)u-^{x))-'"\ Ulix) = {U,^i{x),U,^2{x),U,,3{x)), and U,,M = 
/3+^u~^(a;) — (3^^u~^^{x). To test the Hessian for negative definiteness, we evaluate the quadratic 
form for an arbitrary vector w. We have 

w'Hw = -^c^(a;) (i«iC/^,i(a;) +u;2C/^,2(a;) + W3C/^,3(^c))^ (34) 

z 

which is < for all to, and generically negative unless all terms in the sum vanish. 
Hence, there is a maximum at y, and, together with Eq. (|28|l . we have 

r = sup \r{y{x)) - g{y{x))\ = r(y) - .g(y) . (35) 
In the case ykixk) — x^, this is the maximum principle as stated in Eq. H18|l . 

B Proof of the maximum principle for infinite sequence 
length 

The proof of the maximum prin ciple in the case i V o o closely follows the corresponding proof 
in the two-state model, compare lHermisson et al\ l)2002(l . The idea is to establish upper and lower 
bounds for a system with finite A'^, which can be shown to converge towards each other in the 
limit as — > oo. 

In order to obtain a lower bound, we look at the system locally. To be specific, we consider 
a volume 14, do in the mutational distance space around do, containing the sequences that can be 
reached from sequence do in at most s mutational steps. If this volume intersects a region where 
r has a jump, we take as Vg^da only that part containing do where r is continuous. The number 
of sequences contained in Vs^do is denoted by ^.{Vs^do)- 

The part of the time-evolution operator associated with this volume consists of those ma- 
trix elements where both the row and the column index correspond to sequences in Vg^do- This 
n(Vs,do) ^ 'T.(V^s,do)"dimensional submatrix describes a system with an effective mutational outflow, 
and hence the local growth rate is lower than the global growth rate. Thus, the largest eigenvalue 
Tg^do of this submatrix yields a lower bound for the largest eigenvalue rjv of the whole (but still 
finite) system. 

In order to obtain estimates for rs,do: we use the symmetrised system described by H. This 
can be done because the eigenvalues of the corresponding submatrices are the same. 

We evaluate Rayleigh's principle for the quadratic form for the vector y* = (1,1,...,1) and get 
as lower bound for the largest eigenvalue of the whole system 

_ . _ y^Hy Y.ij{Hs,do)rj 

TN > rs,do = sup — - — > — r — . (36) 

y yty n{V,,do) 

To write this more explicitly, but in a compact way, we introduce the function 



gN.d - J2 + - V "^-^''''^ - y^die.^d^ ) , (37) 



which sums over all mutational terms in the row labelled by d of the full matrix H. Using this, 
we get as a lower bound 

rN > rs,do > —777 — T I ("^d.- 9N.d) + (boundary terms) 1 . (38) 
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The boundary terms are the terms describing the mutational flow through the surface Sg^do of the 
volume Vs^do i^^to Vg^do from the outside, which are contained in the full H, but not in H^do- 



They are of the form \J u^_^^u^^ with d e Vg^do and d + ^ ^s,do- 

To perform the limit ^ cxi as described in section lF31 let ~ r[xd) and = u^^[xd) be 
given as continuous functions, and analogously, gN.d — 9N{xd)- The size of Vs^da shall be scaled 
such that sn ^ \/N. With increasing N, the mutational distances x of neighbouring sequences 
approach each other, and so do the values of the functions r, u^^ and as they are continuous 
functions. More precisely, the total Hamming distance X]fe=i — d'k \ between any two sequences 
d, d' £ Vs.df, is at most 2s. For the mutational distances we then have Xd — Xd' = > with 

increasing N because J2^^ < 2s _ ^-1/2 _^ q. For every x in the mutational distance 

space, we can choose a suitable sequence (dpf) = (d]^(x)) such that Xd,^ x. For any distance d , 
such that dN - d' lies in K^^dw, we then have limAr^oo [r{xd^+d') - gN{xdM+d')] = r{x) - g{x). 

On the other hand, the number of sequences in the volume Vsj^^do increases with N^^'^ , whereas 
the number of sequences in the surface Ss^^do of the volume, and likewise the number of surface 
terms in Eq. (|39|l . only increases with N. Therefore, we get \imN^oc^sf,,dN > r{x) — g{x) for 
arbitrary x. 

For an upper bound, we consider a global maximum of the ancestral distribution, i.e., a 
such that > ad for all d. An evaluation of Eq. (|23|l for d^ and y/od/ < y^ad+ yields 

TN < - gjv,d+ < sup (r-d - gN.d) ■ (39) 

d 

Performing the limit in the same way as above, we get r^+—gj^ ^+ —>■ r{x^)—g(x^), and combining 
this with the lower bound, we have 

sup [r{x) - g{x)] < Too < r{x+) - g{x+) < sup [r{x) - g{x)] , (40) 

X X 

which proves Eq. lfTH|l . 

Now, it remains to be shown that the ancestral distribution is peaked around a;+, and that 
indeed x^ = x as well as r = r(x). For this, we start again with the eigenvalue equation in 
ancestral form 123|l . multiply by -^0^ and sum over the mutational distance space 



rN 



E 



d L 4 



(41) 



Using y^ad±e^ad < \ (ad±e; + Od) , we get 

TN < {rd - gN.d) ad ■ (42) 

d 

As rN — *■ Too and gN.d — * g{xd) uniformly, for every e > we can find an A^^ such that for every 
N>N, 

Too - < ^ {r{xd) - g{xd)) ad ■ (43) 

d 

Due to Eq. H40I) . we have r{xd) ~ g{xd) < ^oo- Splitting the sum into two parts, X]d> +Sd< with 

rixd-^) - g{xd^) >roo-e and r{xd^) ~ gixd<) < t'oo - e, yields 

roo-e^ < ^00 ^ ad:, + (^00 - e) ^ ad< = ^00 - e ^ ad< • (44) 

d> d< <i< 
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Thus, we have X]d< '^d< < e, which means that for N sufRciently large, only sequences with 
r{x) — g{x) arbitrarily close to its maximum contribute to ancestral means. Thus, in the 
generic case that the maximum is unique, the ancestral distribution is peaked around a;+ = x, 
and thus r — r{x), which implies Eq. H20|l . 
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